查看原文
其他

coco 2018-06-06


原来你是我最想留住的幸运,遇到学生信的你们就是最大的幸运,感谢帮助,感谢照顾。


今天看我的基因组27时,猛的发现,我忘了vcf格式了。之前看过的。忍不住反思自己是不是金鱼,7秒的记忆!!


vcf是用于描述SNP,INDEL和SV的结果的文本文件。

samtools和GATK得到的结果文件有一点差别。


CHROM:       表示变异位点是在哪个contig 里call出来的,如果是人类全基因组的话那就是chr1…chr22,chrX,Y,M。

POS:          变异位点相对于参考基因组所在的位置,如果是indel,就是第一个碱基所在的位置。

ID:            如果call出来的SNP存在于dbSNP数据库里,就会显示相应的dbSNP里的rs编号。

REFREF:     在这个变异位点处,参考基因组中所对应的碱基和研究对象基因组中所对应的碱基。

QUAL:         可以理解为所call出来的变异位点的质量值。Q=-10lgP,Q表示质量值;P表示这个位点发生错误的概率。因此,如果想把错误率从控制在90%以上,P的阈值就是1/10,那lg(1/10)=-1,Q=(-10)*(-1)=10。同理,当Q=20时,错误率就控制在了0.01。

FILTER:       理想情况下,QUAL这个值应该是用所有的错误模型算出来的,这个值就可以代表正确的变异位点了,但是事实是做不到的。因此,还需要对原始变异位点做进一步的过滤。无论你用什么方法对变异位点进行过滤,过滤完了之后,在FILTER一栏都会留下过滤记录,如果是通过了过滤标准,那么这些通过标准的好的变异位点的FILTER一栏就会注释一个PASS,如果没有通过过滤,就会在FILTER这一栏提示除了PASS的其他信息。如果这一栏是一个“.”的话,就说明没有进行过任何过滤。

#CHROM  POS ID      REF ALT QUAL    FILTER  INFO    FORMAT  NA12878 chr1    873762 . T G 5231.78 PASS AC=1;AF=0.50;AN=2;DP=315;Dels=0.00;HRun=2;HaplotypeScore=15.11;MQ=91.05;MQ0=15;QD=16.61;SB=-1533.02;VQSLOD=-1.5473 GT:AD:DP:GQ:PL 0/1:173,141:282:99:255,0,255

chr1    877664 rs3828047 A G 3931.66 PASS AC=2;AF=1.00;AN=2;DB;DP=105;Dels=0.00;HRun=1;HaplotypeScore=1.59;MQ=92.52;MQ0=4;QD=37.44;SB=-1152.13;VQSLOD= 0.1185 GT:AD:DP:GQ:PL 1/1:0,105:94:99:255,255,0

chr1    899282 rs28548431 C T 71.77 PASS AC=1;AF=0.50;AN=2;DB;DP=4;Dels=0.00;HRun=0;HaplotypeScore=0.00;MQ=99.00;MQ0=0;QD=17.94;SB=-46.55;VQSLOD=-1.9148 GT:AD:DP:GQ:PL 0/1:1,3:4:25.92:103,0,26

chr1    974165 rs9442391 T C 29.84 LowQual AC=1;AF=0.50;AN=2;DB;DP=18;Dels=0.00;HRun=1;HaplotypeScore=0.16;MQ=95.26;MQ0=0;QD=1.66;SB=-0.98 GT:AD:DP:GQ:PL 0/1:14,4:14:60.91:61,0,255

来来来,说一说栗子:

chr1:873762 是一个新发现的T/G变异,并且有很高的可信度(qual=5231.78)。

chr1:877664 是一个已知的变异为A/G 的SNP位点,名字rs3828047,并且具有很高的可信度(qual=3931.66)。

chr1:899282 是一个已知的变异为C/T的SNP位点,名字rs28548431,但可信度较低(qual=71.77)。

chr1:974165 是一个已知的变异为T/C的SNP位点,名字rs9442391,但是这个位点的质量值很低,被标 成了“LowQual”,在后续分析中可以被过滤掉。

 一看这么复杂


但是里面打都是tag,记不住就查吧(我也很无奈)。其实最关键的信息也就是那么几列:

chr1    873762      .       T   G   [CLIPPED]  GT:AD:DP:GQ:PL    0/1:173,141:282:99:255,0,255

chr1    877664  rs3828047   A   G   [CLIPPED]  GT:AD:DP:GQ:PL    1/1:0,105:94:99:255,255,0

chr1    899282  rs28548431  C   T   [CLIPPED]  GT:AD:DP:GQ:PL    0/1:1,3:4:25.92:103,0,26


其中最后面两列是相对应的,每一个tag对应一个或者一组值,如:

chr1:873762,GT对应0/1;AD对应173,141;DP对应282;GQ对应99;PL对应255,0,255。


GT:    表示这个样本的基因型,对于一个二倍体生物,GT值表示的是这个样本在这个位点所携带的两个等位基因。0表示跟REF一样;1表示表示跟ALT一样;2表示第二个ALT。当只有一个ALT 等位基因的时候,0/0表示纯和且跟REF一致;0/1表示杂合,两个allele一个是ALT一个是REF;1/1表示纯和且都为ALT; The most common format subfield is GT (genotype) data. If the GT subfield is present, it must be the first subfield. In the sample data, genotype alleles are numeric: the REF allele is 0, the first ALT allele is 1, and so on. The allele separator is '/' for unphased genotypes and '|' for phased genotypes.

0 - reference call

1 - alternative call 1

2 - alternative call 2

AD:    对应两个以逗号隔开的值,这两个值分别表示覆盖到REF和ALT碱基的reads数,相当于支持REF和支持ALT的测序深度。

DP:    覆盖到这个位点的总的reads数量,相当于这个位点的深度(并不是多有的reads数量,而是大概一定质量值要求的reads数)。

PL:      对应3个以逗号隔开的值,这三个值分别表示该位点基因型是0/0,0/1,1/1的没经过先验的标准化Phred-scaled似然值(L)。如果转换成支持该基因型概率(P)的话,由于L=-10lgP,那么P=10^(-L/10),因此,当L值为0时,P=10^0=1。因此,这个值越小,支持概率就越大,也就是说是这个基因型的可能性越大。

GQ:   表示最可能的基因型的质量值。表示的意义同QUAL。



举个例子说明一下:


chr1    899282  rs28548431  C   T   [CLIPPED]  GT:AD:DP:GQ:PL    0/1:1,3:4:25.92:103,0,26

在这个位点,GT=0/1,也就是说这个位点的基因型是C/T;GQ=25.92,质量值并不算太高,可能是因为cover到这个位点的reads数太少,DP=4,也就是说只有4条reads支持这个地方的变异;AD=1,3,也就是说支持REF的read有一条,支持ALT的有3条;在PL里,这个位点基因型的不确定性就表现的更突出了,0/1的PL值为0,虽然支持0/1的概率很高;但是1/1的PL值只有26,也就是说还有10^(-2.6)=0.25%的可能性是1/1;但几乎不可能是0/0,因为支持0/0的概率只有10^(-10.3)=5*10-11


下面介绍一下第8列信息:




该列信息最多了,都是以 “TAG=Value”,并使用”;”分隔的形式。其中很多的注释信息在VCF文件的头部注释中给出。以下是这些TAG的解释

AC,AF 和 AN:AC(Allele Count) 表示该Allele的数目;AF(Allele Frequency) 表示Allele的频率; AN(Allele Number) 表示Allele的总数目。对于1个diploid sample而言:则基因型 0/1 表示sample为杂合子,Allele数为1(双倍体的sample在该位点只有1个等位基因发生了突变),Allele的频率为0.5(双倍体的 sample在该位点只有50%的等位基因发生了突变),总的Allele为2; 基因型 1/1 则表示sample为纯合的,Allele数为2,Allele的频率为1,总的Allele为2。

DP: reads覆盖度。是一些reads被过滤掉后的覆盖度。

Dels: Fraction of Reads Containing Spanning Deletions。进行SNP和INDEL calling的结果中,有该TAG并且值为0表示该位点为SNP,没有则为INDEL。

FS:使用Fisher’s精确检验来检测strand bias而得到的Fhred格式的p值。该值越小越好。一般进行filter的时候,可以设置 FS < 10~20。

HaplotypeScore: Consistency of the site with at most two segregating haplotypes

InbreedingCoeff: Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hard-Weinberg expectation

MLEAC: Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), for each ALT allele, in the same order as listed

MLEAF: Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT alle in the same order as listed

MQ: RMS Mapping Quality

MQ0: Total Mapping Quality Zero Reads

MQRankSum: Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities

QD: Variant Confidence/Quality by Depth

RPA: Number of times tandem repeat unit is repeated, for each allele (including reference)

RU: Tandem repeat unit (bases)

ReadPosRankSum: Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias

STR: Variant is a short tandem repeat 

点击阅读原文,有

我的基因组28-必须要理解vcf格式记录的变异位点信息

还有更多文章,请移步公众号阅读

如果你生信基本技能已经入门,需要提高自己,请关注上面的生信技能树,看我们是如何完善生信技能,成为一个生信全栈工程师。

如果你是初学者,请关注下面的生信菜鸟团,了解生信基础名词,概念,扎实的打好基础,争取早日入门。

    您可能也对以下帖子感兴趣

    文章有问题?点此查看未经处理的缓存